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(ABSTRACT) 


Numerical solutions are obtained for the quasi-compressible Navier- 
Stokes equations governing the time dependent natural convection flow 
within a horizontal cylinder. The early time flow development and wall 
heat transfer is obtained after imposing a uniformly cold wall boundary 
condition on the cylinder. Solutions are also obtained for the case of 
a time varying cold wall boundary condition. Windward explicit differ- 
encing is used for the numerical solutions . The viscous truncation 
error associated with this scheme is controlled so that first order 
accuracy is maintained in time and space. The results encompass a 
range of Grashof numbers from 8.34 x 10 4 to 7 x 10 which is within the 
laminar flow regime for gravitationally driven fluid flows. Experiments 
within a small scale instrumented horizontal cylinder revealed the time 
development of the temperature distribution across the boundary layer 
and also the decay of wall heat transfer with time. Agreement between 

measured temperature distributions and the numerical solutions was 

l/k 

generally good. The time decay of the dimensionless ratio Nu/G r is 
found numerically and experimentally and, over most of the cylinder wall, 
good agreement is obtained between these two results . The numerical 
results indicate that the fluid exhibits a strong tendency to resist 
first order motion within the inner core region. The early establish- 
ment of a shallow positive upward temperature gradient within the core 
enhances its stability. No first order vortical motion is induced by 
the boundary layer and this is attributed in part to the fluid decelera- 
tion near the bottom of the cylinder along with expulsion of fluid from 



the boundary layer in the lower portions of the cylinder. 
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INTRODUCTION 


Chapter I 

Natural convection flows within closed containers have intrigued 
mathematicians and fluid dynami cists for many years. Part of the 
motivation for understanding such flows was based upon an early 
realization that many practical engineering situations are governed by 
gravitationally driven fluid flows. Present day interest is concerned 
with the flow of fluids within pipes , nuclear reactor cooling systems , 
turbine blades, and stationary containers. These fluid flows may be 
significantly influenced by natural or induced body forces. The 
resultant fluid motion and heat transfer are the principle features of 
interest and it is toward an understanding of these features that most 
studies have been directed. The early work by Nusselt (reference l) , 
Hermann (reference 2), Beckmann (reference 3), and Hermann (reference 
4) established both the appropriate governing equations as well as the 
general relationship between the non-dimensional heat transfer and the 
Grashof number for external flows over cylindrical configurations. 
Numerous experimental studies have confirmed the theoretical findings 
(reference 5 through 7) for large values of the Grashof number wherein 
a boundary-layer flow is present. Not until the work reported by 
Ostrach (reference 8), Lewis (reference 9), Batchelor (reference 10), 
and Pillow (reference 11 ) was the internal natural convection problem 
given a general theoretical treatment comparable to the external 
problem. More recent analytical studies involving somewhat restricted 
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wall boundary conditions have been reported by De Vahl Davis 
(reference 12), Weiribaum (reference 13), Menold (reference lU), Hantman 
and Ostrach, (reference 15) and Gill (reference l6). 

The principal problem that has thus far prevented a general 
analytical solution stems from the coupling between the boundary 
layer flow near the walls of the container and the core flow that is 
driven by the boundary layer flow. In formulating the problem in terms 
of a stream function most of the previous investigators except Hantman 
and Ostrach (reference 15) have been faced with the necessity of 
specifying a core stream function behavior that is not known a priori. 
Thus , depending upon specified wall boundary conditions , the core 
has been assumed to have either an isothermal constant vorticity 
character or to be stratified with streamlines extending into the 
boundary layer. Both analytical solutions and experiments have shown 
that two entirely different flow configurations are possible. If the 
flow is heated from below, the core streamlines are closed and the core 
is isothermal. If the wall boundary condition is such that heating 
occurs from the side, the core is stratified with isotherms and 
streamlines coinciding. The studies reported thus far are not able to 
predict the critical heating phase angle at which a rotating flow 
configuration changes to a stratified flow configuration. The Oseen 
linearization used by Weinbaum does not help the problem because this 
approximation decouples the core flow from the boundary layer flow. 

The core is expected to be closely coupled to the boundary layer flow 
because it is driven by the boundary layer. 
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A related problem with less troublesome boundary conditions is 
that of a horizontal cylinder with one half of its bounding walls 
raised to a constant, uniform temperature and the other half lowered 
an equal amount below the initial fluid temperature. Thus a rotating 
boundary-layer flow encircles the inside cylindrical walls with the 
flow rising over one half the circumference and falling over the other 
half. This problem was studied numerically by Heliums (reference 16). 
Heliums used what is now termed the windward or donor cell differencing 
technique for the time dependent problem starting from initially imposed 
wall boundary conditions. He obtained numerical solutions for the 
velocity, temperature, and heat transfer distributions within the 
cylinder. Because fluid is both entrained and ejected by the boundary- 
layer flow, the windward finite differencing scheme is quite suitable 
for this type of problem. Heliums was able to make favorable compari- 
sons between his solutions and experiments carried out by Martini and 
Churchill (reference 17). These comparisons were possible for the 
steady state flow only. No unsteady flow measurements were reported. 

For steady state flow, Heliums verified the relationship between the 

lA 

Nusselt number and Grashof number; Nu » C G ' both from a formal 
derivation of the dimensionless governing equations and from the 
resultant numerical solutions. For unsteady flow the coefficient C 
can be expected to be time dependent. The work reported by Heliums 
is closely related to the present study and will be discussed in more 
detail in the chapters following. At the present time only a very few 
of the possible boundary conditions that could be imposed on a 
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horizontal cylinder have been studied and with only partial success 
in most of the cases reported. 

The work reported here represents an attempt to clarify the flow 
resulting from a new wall boundary condition that has not been 
here-to-fore studied for the internal flow problem. A uniform cold 
wall is established at time zero and the early time development of the 
fluid motion is studied. The fluid is initially at rest and will also 
return to rest at very large times after flow initiation. Thus a 
steady state is not of interest in the present work. In addition, 
solutions are obtained for the case in which the wall temperature 
decays with time. This problem corresponds to experimental boundary 
conditions that were imposed within a small cylinder in the present 
work. The experiments are described and a discussion of the relation- 
ship between the numerical work and the experimental work is given. 



A Physical Description of the Problem 


Chapter II 

The geometry and nomenclature of the horizontal cylinder to be 
studied is shown in figure 1. The cylinder is of semi infinite length 
to allow a region of two dimensional flow to exist^. One practical 
application for such a geometry is related to a large blowdown wind 
tunnel storage facility in which the major heat loss to the walls 
occurs due to gravitationally driven natural convection flows. When a 
hot gas is stored in such a tank appreciable azimuthal gas flow occurs 
due to the imbalance between the gravitational body forces and the 
existing hydrostatic pressure gradient in the fluid. The gas, which is 
air, is initially at rest with a balance between the body force and the 
hydrostatic pressure gradient. The gas is initially at a uniform 
temperature T^,. At time zero, a uniform cold wall is imposed on the 
cylinder, and the resulting conduction of heat out of the gas near the 
wall causes the convective motion to begin . As the flow develops , a 
thin boundary layer is formed near the wall and this layer thickens 
with time. The initially motionless inner fluid (the core) is driven 
by the boundary layer flow and gives up energy to the heat conducting 
boundary layer fluid. As time progresses the boundary layer will 
affect the inner-most regions of the core flow and eventually the gas 
can be expected to give up all of its excess energy to the cold wall. 

1 A two dimensional flow within a horizontal cylinder has been observed 
by Brooks and Ostrach (reference 21 ). 
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When this happens the fluid will he once again at rest with a uniform 
temperature now equal to the wall temperature. Several distinct fea- 
tures are intuitively apparent. Because of the uniform wall tempera- 
ture, a mid-plane of symmetry is immediately established with the 
dividing line running vertically upward through the center of the 
cylinder. Along the mid plane of symmetry the azimuthal velocity is 
zero. In each half cylinder there are two stagnation points which will 
be at the intersection of the line of mid-plane symmetry and the wall. 

The principle motion in the boundary layer will be azimuthal and 
the induced motion will be radial. Because of the small coefficient 
of viscosity for air we may expect velocities of lower order in the 
core then in the boundary layer as well as both inward and outward 
radial velocities across the boundary layer. In addition it might be 
anticipated that some stratifipation of the flow will occur in the 
lower portions of the cylinder. The geometry and principle features 
of the flow having been outlined, the remainder of this thesis will be 
concerned with obtaining an understanding of the details of the flow 
from numerical solutions to the governing equations as well as 
experimental measurements. 



Mathematical Formulation of the Governing Equations 


Chapter III 

The model chosen for this study is that of a viscous, heat 
conducting, quasi-compressible fluid that conforms to the Bous s inesq 
approximation. For small differences between the gas temperature and 
the wall temperature, the density may be taken as a function of 
temperature only and considered as a variable only where it modifies 
the body force terms in the equations describing conservation of 
momentum. This approximation has been investigated extensively in 
references 4, 8, 15 » and l6. 

In dimensional form the quasi -compressible Navier Stokes equations 
applicable to the case of large Grashof Numbers and small gas to wall 
temperature differences in cylindrical coordinates are: 
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Energy 


3.U 


Viscous dissipation and pressure work contributions to the energy 
equation are neglected at the outset as these terms have been shown 
to be negligible for gravitationally driven flows. (See Ostrach 
reference 10). 


P 



BousSinesq Equation of State for 
a perfect gas 3.5 


This system of equations consists of three second order, non- 
linear partial differential equations, one first order, linear partial 
differential equation and an algebraic equation of state. For the 
physical flow described in chapter two, the following initial and 
boundary conditions axe imposed: 


at t = 0 

u(r,0,O) = v(r,0,O) = 0 
T(r,0,O) = 

p(r,0,O) = p^r.e.O) 
p(r,0,O) = P i 


3.6 
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at 


r = 


w 


u(r w ,0,t) = v(r w ,6,t) = 0 
T(r w , 0 ,t) = T w 


and at r = 0 


p(o,e,t) = p q 


and on 


6=0 and v 



Mid Plane Symmetry 3.8 



3T 

36 


0 


The equations 3.1 - 3.8 will he more conveniently dealt with in 
non dimensional form. The following dimensionless quantities are 
defined: 


U = 


u 

/g 8a r w 


V 

8a \ 


V = 
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Substituting equations 3.9 into equations 3-1 - 3.8 and rearranging 
gives : 
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U = V = 0 

<t> ■ 1.0 


3.15 
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and at R - 0 

P = P 

o 

on 0 = 0 and it 

QL - ^ = ~j-- = 0 Mid Plane Symmetry 3.17 

90 90 36 

The nondimensionalization provides the unit order variables U, V, <t> » 
and p. The reference pressure is Pjg Ba r w so that the dimension- 
less pressure will he greater than order one. The reference time is of 
order one, so the dimensionless time can also he greater than unity. 

In previous investigations some simplification of equations 3*10 
through 3.13 were made. Heliums (reference l6) derived the governing 

equations for unsteady internal flows within horizontal cylinders. 

3V 

The Coriolis force as well as the derivative in the azimuthal 

momentum equation were shown to be negligible. In addition, the 

viscous terms in the radial momentum equation and the first three 

terms on the left hand side of equation 3.11 are negligible, (see 

Hermann, reference 1+). Finally, the derivative “ was found to he 

80 

a low order term in equation 3.13. However the derivatives considered 
negligible were in fact carried along within the present computational 
scheme and systematic checks were made to verify that the small order 
of the terms persisted. Formal ordering of equation 3*10 through 3*1^ 
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by asymptotic series expansions of the dependent variables and stretch- 
ing of the normal coordinate is complicated by the eventual thickening 
of the boundary layer. At large values of time the boundary layer is no 
longer thin and the matched core expansions and inner expansions are 
invalidated. An additional simplification to equations 3.10 through 
3,14 involves the form of the pressure gradient terms. 

It is convenient to divide the pressure into two separate parts: 


P * P i + P* 


3.18 


P. is the initial zero motion hydrostatic pressure. P* is a "dynamic" 
pressure that arises due to motion and thermal changes within the 
fluid. The motivation for equation 3.18 is well founded in that the 
initial hydrostatic pressure gradient is known exactly: 


1 3P j _ sin 6 
pR 36 fT~ 


3.19 


1 3P i 1 
i = 4 cos 6 
p 3R p 


3.20 


For many gravitationally driven flows the contribution to the overall 

3P# 3P* / 

pressure gradient from ^g -- or a order effect, (see 

Heliums reference 1 6 and Ostrach reference 18). Knowledge of the 

hydrostatic pressure gradient facilitates solution of the governing 

equations through the use of equations 3*19 snd 3*20. With these 

considerations and substituting equations 3*19 and 3*20 into equations 



lU 


3.10 through 3.1*+ the governing system becomes: 
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3.25 


at 


R = 1.0 
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U = V = 0 3.26 

<t> = 1.0 

and at R = 0 

P* = P* 
o 

and on 0 — 0 and tt 


iSU o 

90 96 


3.27 



Development of the Finite Difference Approximation 
to the Governing Equations 

Chapter IV 

Low velocity fluid flows have proven difficult to solve using 
finite difference techniques. When the fluid velocity is low, long 
computational times are required, and for gravitationally driven flows 
that start from rest, the situation is further aggravated. The fluid 
flow studied here should develop maximum velocities on the order of 
several feet per, second; thus long computational times appear unavoid- 
able. The second difficulty arises from the fact that both positive 
and negative velocity components may be expected as fluid will be both 
entrained and expelled by the boundary layer. Numerical instabilities 
are produced by positive and negative coefficients of the convective 
terms of a difference scheme unless special care is taken regarding the 
form of the differencing. For these reasons an explicit, windward 
differencing technique was chosen (also called the donor cell technique). 
As shown recently by Roache (reference 19), the windward scheme along 
with all but one other scheme (reference 20) have only first-order 
accuracy in the transient development of the flow. A diffusive trunca- 
tion error is introduced which appears as an artificial viscosity that 
may override the physical viscous damping within the fluid. Because of 
this, care must be taken with the windward scheme to insure that the 
grid size is sufficiently small so that artificial viscous effects do 
not have a predominant effect on the numerical development of the flow. 
From reference 19 it is seen that the artificial viscosity of the 

16 
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transient windward differencing scheme in dimensional form is : 


- l ulAr ( , hJ»At v 

~ 2 V Ar' 


4.1 


If a is to be much less than the physical viscosity we may write: 
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lujAr / n l ulAt 
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4.2 


° r to «ai + l ui4t 

lul 

Providing that 

u At « 7^7 4.3 

PM 


Thus , depending on the fluid velocity and the allowable time step 
obtained from stability constraints, the step size must be kept well 
below the right hand side of inequality 4.3. For some flows that have 
been studied, this requirement is over restrictive and can be relaxed. 
(See Callens , reference 22)'. With these considerations the finite 
difference equations approximating eqs. 3.18 - 3.25 are written as: 
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at 


R = 0 


P* = 0 


Because both U and $ are very small quantities near R = 0 
the boundary condition for P* given in equations 4.9 is a close 
approximation to the physical situation. 

The stability and convergence of this system of equations is 
inferred by applying a Von Neumann stability analysis to the 
linearized forms of the difference equations. Appendix A provides the 
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details of the analysis. Although the procedure outlined by Richtmyer 
(reference 26) is rigorously applicable to a limited class of linear 
equations, it has worked successfully for the non-linear equations 
of fluid mechanics. The allowable time step obtained in appendix A is 
given by: 


At < 


- 2 


R(Jt)A6 


AR y$T (R(p)A6) 2 AR £ 


4.10 


The time step chosen for most of the calculations was 0.8 of the 
right hand side of equation 4.10. The value provided stable calcula- 
tions over the majority of cases studied. Theoretically it is difficult 
to establish optimum grid spacing ratios. Past work across boundary 
layers has been dictated by the need to resolve large gradients over 
short spatial lengths, and thus a fine grid spacing has been used in 
the spatial dimension normal to a bounding surface. For the cylindrical 
geometry the principal flow in the boundary layer is azimuthal, but it 
was found that radial flow is of primary importance in the core. A grid 
network that is refined in only the radial direction (normal to the 
cylinder walls) does not seem appropriate here. The computational 
experience with different grid spacing ratios indicates that a ratio 

— << 1.0 is essential for stable results. This probably should not be 
A0 

surprising in a cylindrical coordinate system because the azimuthal 
grid expression appears as R(JI)A0 in the difference equations, and it 
is immediately seen that this value reduces to ARA6 at the inner most 
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grid location next to the origin of coordinates. 

The for mal accuracy of the windward difference scheme is related 
to the computational cell Reynolds number. For the present study this 
number may be written as: 


R . = IUI ar/T U.ll 

eA r 

For R « 2 the windward differencing will have first order accuracy 
ea 

in both the spatial and temporal differences. 



V Numerical Solutions 


A. Constant Wall Temperature 

The system of equations 4.4 - 4.9 was solved for values of U, <j», 
V, and P* over the entire grid following a step- function change in 
the cylinder wall temperature at time zero. A flow diagram describing 
the calculation sequence is shown in figure 3. It is seen that the 
correct difference equations are selected at each grid point depending 
upon the sign of U. . or V. - such that stable computations will 
result. The program was set up so that computations could be continued 
from a previous computer run, and thus extended run times up to 5 hours 
on the Langley Research Center CDC 6600 computer were made possible. 
Typical computational speeds for a 51 by 31 grid for the half cylinder 
were 1.02 x 10^ grid points per hour. The real time development of 
the flow progressed at a rate of 1.25 seconds per hour of machine time. 
An order of magnitude faster flow development time was obtained with a 
26 x 16 grid network; however a formal first order accuracy is not 

achieved with such a coarse grid. 

The steps taken in the computational sequence of events were as 

follows : 

1. Equation 4.7 was solved over the entire grid for new values 
of <{i. 

2. Equation 4.4 was then solved for new values of U over the 
entire grid. 

3. Values of <j> and U were updated. 
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4. Equation 4.6 was solved for the current values of V over the 
entire grid. 

5. Equation 4.5 was solved for current values of P* over the 
entire grid. 

6. Values of V and P* were updated. 

7. Equation 4.10 was solved to determine an allowable time step. 

8. Values of the Nusselt number at the wall were calculated. 

The process was repeated over the entire grid for the next time step. 
The dimensionless heat transfer at the wall is given by the Nusselt 
number and obtained from the following: 


or 


so that 
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or in dimensionless form: 


s (1 - V waU -i> 


5 . 1 * 


Table I lists the values of the input parameters that were used for 
the solutions to be presented. Figure 4 shows the results obtained for 
the case II solutions of Table I. 
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From the simple physical geometry of the cylinder a rather compli- 
cated flow pattern is revealed by the numerical results. Figure 4a shows 
the azimuthal velocity distribution at 6 - ir/2 for three different 
time levels. The early time T = 3.05 can be considered as showing the 
distribution prior to what is considered to be fully developed flow. 

By this it is meant that peak velocities at any given location have not 
been achieved. The intermediate time, T = 4.57, shows the velocity 
distribution at its peak value, and the late time, T = l4. 5, shows the 
distribution after a decay in the velocity has taken place. Significant 
inward displacement of the peak velocities is not apparent. The 
distributions do thicken with time over about 20 grid point spacings, 
and the viscous effects are transported further into the core of the 
fluid as time progresses • The distributions shown in figure 4 are 
typical of all the results obtained for the constant wall temperature 
case. Figure 4b shows the velocity distribution near the bottom of 
the cylinder at an azimuthal angle of 23.7°. The momentum gathered by 
the fluid falling downward has both thickened the boundary layer as 
well as increased the peak azimuthal velocity. At values of 9 less 
than 23° the fluid rapidly decelerates and comes nearly to rest in 
the lower part of the cylinder. Figure 4c shows the velocity distribu- 
tion near the top of the cylinder. The boundary layer is well 
defined, but the azimuthal velocities are substantially lower than in 
the bottom portion of the cylinder near the wall. Figure 5a shows the 
radial velocity distribution for case II of Table I. Here it is seen 
that the peak velocities are an order of magnitude less than 
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corresponding azimuthal velocities. The radial velocities peak at 
their maximum values at the later time, T = 1^.5, in contrast to the 
azimuthal velocity peaks which occur near T = 4.57* At 0 = 23.7°, 
fi gixre 5b, the radial velocity is negative which indicates an explusion 
of fluid from the boundary layer. As the lower portion of the cylinder 
"fills" with fluid that has moved downward in the boundary layer, 
negative or radially inward flow takes place to accommodate the added 
fluid. Thus, fluid is slowly forced into the boundary layer near the 
top of the cylinder due to displacement effects from below. Figure 5c 
shows the radial velocity distribution in the upper portion of the 
cylinder. At 6= l6l° the flow is still developing with peak 

velocities occurring near T = lU . 5 • 

Thus the general picture of the flow field involves three major 
features . First, there is a boundary-layer development near the wall 
due primarily to downward azimuthal fluid flow near the wall. Secondly, 
an induced radial velocity occurs that feeds fluid into the developing 
boundary layer in the upper and middle azimuthal locations of the 
cylinder and ejects fluid out of the boundary-layer at lower azimuthal 
locations. Third, a core region exists that strongly resists first- 
order motions and only very slowly forces fluid at lower levels to 
rise and enter into the boundary-layer . This is a striking example 
of a flow in which the principal motion is confined to the boundary 
layer . 

The dimensionless temperature, T/T^, is shown as a function of 
radial distance in figures 6. Here the thermal-boundary-layer is less 
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well developed than the velocity-layer at a given time. The conduction 
terms in the energy equation cause this behavior due to the very low 
fluid velocities that occur at early times after the cold wall initial 
condition is imposed. The accuracy of the results increases as the grid 
network is refined, but the computational time required for fine grids 
is enormous. The early behavior of wall— heat transfer is of principle 
interest, and figure 7 shows the Nusselt number decay with time for 
three different azimuthal stations. The computational results show that 
heat is conducted out of the fluid at early times at a greater rate than 
energy is convected into a fluid element. As the velocity field 
develops, this trend is altered causing a slight steepening of the 
temperature gradient near the wall. This result appears to be a valid 
physical description of the flow development. The convective terms in 
the energy equation are negligible when very low fluid velocities are 
present during early times. The result is an energy balance that is 
dominated by conduction out of the fluid to the wall until the flow 
field develops. 

The increasing value of Nusselt number with increasing 6 gives a 

clear picture of the positive upward temperature gradient within the 

boundary layer. This positive, upward gradient also exists in the 

core fluid as indicated by figure 8. The early establishment of this 

upward gradient produces a thermally stable core which tends to resist 

downward motion. The low velocities calculated in the core flow are in 

N 

part, a result of this thermal behavior. Figure 9 shows a plot of U 
as a function of time. The cold fluid entering and residing in the r 


IF* 
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lower portion of the cylinder causes a rapid and nearly linear decay 
of the dimensionless heat transfer function. For 0 near tt/ 2 and up 
to 0 = tt, the function decays less rapidly due to the flow of warm 
fluid into the houndary layer from the core. At very large times the 
heat transfer function must asymptotically approach zero because the 
fluid gives up all of its excess energy to the wall. It is apparent 
that C 1 -*■ 00 as t -► 0 due to the step function cold wall initial 
condition. The actual value of at early times near t - o is 

dependent on the radial grid spacing used, as might be expected. A 
comparison of the effects of grid refinement is shown in figure 10. 

Here the value of T/T^ as a function of radial distance is shown for 
three different grid networks. The radial coordinate grid spacings are 
.02; .01; and .0067. The azimuthal location is 0 = 90° , and the time 
is T = 2.7. 

The profiles steepen as the grid is refined and convergence of 

the finite difference solution is evident. Of interest is the velocity 

distribution shown in figure 10b. Here convergence is also indicated 

by the refined grid results. The coarse grid velocity peak is, 

however, above the finer grid peaks. The slope of the distributions 

near the wall appears to be adjusting itself from an "overshoot" where 

AR = .01 to a convergent result as AR becomes smaller. Similar 

h 7 

results have been computed at Grashof numbers of 8 x 10 and 7 x 10 
for the grid range given above. Case I of table I was computed for a 
Grashof number of 8.3 x lo\ The computational time increases 
considerably at this relatively low Grashof number. Figure 11-a shows 



28 


some representative azimuthal velocity distributions for this case. 

The lower Grashof number results for this case can be interpreted as 
being due to a more viscous fluid flow, and the distributions reflect 
this. At a dimensionless time of T — 3.07 the boundary layer is 
considerably thicker than for the comparable boundary layer thickness of 
case II. Even at the early time of t = 1.31, the viscous effects are 
more pronounced than for the higher Grashof number case . Figure 11-b 
shows typical radial velocity distributions for the low Grashof number 
case. At t — 3.07 the distribution shows the effects of higher 
viscosity with a thicker profile than for the case II solution. Finally, 
figure 11-c shows the dimensionless temperature function distribution. 
Again the more viscous fluid of case I shows a thicker thermal 
boundary layer than the case II solution. All of the profiles for case 
I have a qualitative resemblance to the results for case II. The lower 
Grashof number solutions for the Nusselt-Grashof relation are shown in 
figure 11-d. These curves should be compared with the results shown 
in figure 9 - 

The computations for the dimensionless dynamic pressure, P*, 
indicate that the dynamic pressure gradient makes only a small order 
contribution to the momentum balance within the fluid. This occurrence 
ia consistent with the results reported by Ostrach in reference 18 and 
2h. 

The results of the computations for case III at a Grashof number 
of 1.3 x 10 6 are shown in figure 12. The velocity profiles and 
temperature profiles represent a result that is intermediate to the case 
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I and case II solutions . 


B. Time Varying Wall Temperature 

Comparison of the experiments to be described in Chapter VI with 
numerical solutions necessitated the use of a time-dependent wall 
temperature boundary condition. The experimental wall temperature 
decay can be described by solution to the one-dimensional heat conduc- 
tion equation. . For a wall of thickness x q initially at a temperature 
0(x,o ) = 1 where: 
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And t = 0 


T = T. 
w 1 


$(x,o) = 1 


T x = Temperature of surface at x = x q 
0 after the cold wall boundary 
condition has been applied. 


the temperature history at the inner face x ■ o is given by: 
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With the solution to equation 5*7 we may calculate the time dependence 
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of the coefficient (G ) _1 ^ 2 that appears in equations 4.4, 4.7, and 
4.10 thus: 





1 



(l - $(o,t)) 


5.8 


With this value calculated at each time step, the system of difference 
equations 4.4 through 4.10 can he solved using each new value of Grashof 
number which results from the time-varying wall temperature. As an 
example of the complementary error function. solution for eq. 5*7, 
figure 13 shows the time— dependent dimensionless wall temperature decay. 
Several different limits on the number of terms taken in the summation 
of equation 5.7 showed rapid convergence for values of n greater than 
about 6. In figure 13 n = 9 was used for the temperature function 
solution. 

A fundamental difficulty is encountered using the error function 
solution given by equation 5*7. At very early times the function is so 
close to unity that a near singularity is introduced into equation 5*8. 
The difference scheme is limited, however , to s m a l l time steps for 
solution of the momentum and energy equations. The wall temperature 
decay obtained from equation 5.7 is numerically incompatible with the 
difference scheme. This difficulty can be overcome by approximating 
equation 5.7 at early times by linear functions of the form: 


$(o,t) = a - bx 


5.9 


y 
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An approximation to the error function solution can be made using 
equation 5.9 if we choose: a <_ 0.99 and b = 0.050. Figure lU shows 

the azimuthal velocity distribution for the case where the wall 
temperature decays according to equation 5*9 rearranged: 

T = *(o,t)(T. - T ) + T 5.10 

v 1 X o X o 

The results shown in this figure correspond to the case II input values 
of Table I. 

The time- dependent wall temperature solutions exhibit several 
interesting differences from the constant wall temperature solutions. 

In effect , a slowly cooled wall-boundary condition produces an early 
time driving force that is quite small. This is characterized by the 
time development of the Grashof number. Figure 15 shows a typical time 
history of the Grashof number for the case II input values and a wall 
that follows a temperature history given by equations 5*9 and 5*10. 

The development of the flow field is directly related to the time 
dependent Grashof number. The immediate consequence of a slowly cooled 
wall should appear as a less fully developed flow at any given time than 
one for which a step function wall temperature change has been imposed. 
The early time behavior is closely related to a lower Grashof number 
flow field. The boundary layer development follows the rising Grashof 
number with a different behavior than for the constant wall temperature 
results. In effect the driving force increases with time, and thus 
the flow experiences a longer and more pronounced acceleration. Figure 
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l6-a shows a typical temperature distribution for the time-dependent 
wall boundary condition as well as for the constant wall temperature 
case. The details revealed by these curves are physically valid in 
that the constant wall temperature case, of necessity, must have a 
steeper slope near the wall and, in addition, must have a thinner 
thermal boundary layer due to the higher Grashof number at equal 
values of time. These features are clearly evident in figure l6-a. 

A comparison of the constant and variable wall temperature cases 
at equal values of the Grashof number reveals the fundamental differences 
between the two cases. Figure l6-b shows that the variable wall 
temperature still has the effect of producing a more viscous flow 
by maintaining lower Grashof numbers throughout the entire time 
development than the constant cold wall case maintains. A comparison 
with the constant Grashof number case is difficult because of the 
fundamental difference between the two case histories . If the two 
boundary layer thicknesses are compared at equal values of R and 
time, we find the results shown in figure 17-a. The difference in the 
peak values of the dimensional velocity are of course even greater than 
shown in this figure because of the different value of reference velocity 
used in the non-dimensionalization, i.e. u =/g $cx r^ U. A 
comparison of the dimensionless azimuthal velocity distributions at 
equal values of both Grashof number and time is shown in figure 17-b. 

As expected the time varying wall case shows a thicker profile but 
with lower peak velocities. Similar results were observed for the 
temperature distribution which is also an indirect result of the 
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viscous behavior of the fluid. Perhaps the most significant single 

consideration regarding the energy transfer to the walls is the time 

dependent Nusselt-Grashof correlation. Heliums (reference 16) 

lA 

established the constant C in the relation H u = C . For 

steady flow within a horizontal cylinder C = 0.326. If a frilly 

developed flow occurs in the present problem, the time dependent C(t) 

should approach the steady value prior to a decay in the fluid motion. 

As time increases indefinitely in the present unsteady flow, C (t) -*■ 

o, since the fluid will then have given up all its excess energy to the 

cold walls. Figure 18 shows the time- dependent behavior for the Nusselt- 

Grashof correlation. The relation N is seen to rapidly approach 

u r 

the steady state or fully developed value. The unsteady flow value of 
this relation for the present problem must be asymptotic with zero at 
large values of time. This behavior is demonstrated by the numerical 
results represented in figure 18. 

Figures 19-a through 19-d illustrate the time history of "fluid 
particles" within the flow field from a sequence of photographs which 
were taken from an oscilloscope display of the particle displacements as 
computed from the numerical solutions. As shown in these figures par- 
ticles were located, at t * 0, along rays separated by an azimuthal dis- 
tance of tt/ 8. From R = 0 to R = .64 the particles were spaced radially 
with AR = .04. From R = .64 to R = 1.0 the radial spacing was AR = .01. 

A two dimensional interpolation was used to obtain particle displacements 
between grid locations for the 31 by 101 mesh used for this case. The 
time dependent wall temperature decay for case II-T of table I was used. 
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The sequence of photographs shows the fluid deceleration and strati- 
fication near the "bottom of the cylinder. Regions where the number of 
particles increase represent regions of higher fluid density as might 
he expected. Figure 19-e is a plot of the displacements at 6.6 
seconds. In this figure the particles belonging to each original ray 
have been faired in to illustrate the displacement profiles. Real 
time movies were made from plots such as these to give a physical 
picture of the flow field development. Because of the low velocities 
within the core flow, large values of time must be obtained with the 
numerical solutions before substantial particle displacements can be 


observed within the core. 



Chapter VI - EXPERIMENTAL STUDIES 


Apparatus 

The instrumented stainless steel cylinder used for measuring 
temperature distributions in a gravitationally driven flow field is 
shown in figure 20-a. Experiments were made with the cylinder in a 
horizontal position, and rotation of the cylinder allowed measurements 
to be made at different azimuthal locations. Gage marks and internal 
thermocouple probes were used as references in setting the cylinder 
in a desired position, and a positive lock cradle was made to insure 
that the cylinder maintained a given position. An optical transit was 
used for precise orientation of the entire system. Figure 20-b is a 
schematic diagram of the entire apparatus. The cylinder was made of 
schedule UO non magnetic stainless pipe with 300 series stainless steel 
end caps welded in place. The inside diameter was nominally 5-9^ 
inches and the length was 60 inches. The cooling manifolds were 
insulated from the outer cooling jacket by means of one inch thick 
micarta rings. The purpose of this insulation was to allow filling of 
the coolant tanks prior to a run, and cooling the liquid to a uniform 
temper at tire without introducing conduction effects downward to the 
actual test cylinder. A threaded connection between the outer jacket 
and the test cylinder established a conduction path that had to be 
accounted for, and thus the insulator rings were installed. 

The inside surface of the test cylinder was machined to a smooth 
finish with average surface projections not greater than 220 micro 
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inches from peak to valley. The cylinder vas cleaned and vacuum leak 
tested prior to installation of the instrumentation. The cylinder 
maintained a pressure of less than 10 ^ mm of mercury for a 2k hour 
period. 

Figure 21 is a photograph of the cooling tanks , inlet manifolds , 
and outer' cooling Jacket surrounding the test cylinder. The valves 
for dumping the cooling fluid into the outer cooling Jacket were 
manually operated, and the liquid coolant in the tanks could he . 
completely discharged in less than one second. The "o" ring sealed 
valves shown in figure 20-h proved exceptionally reliable, and because 
of their large size the effective discharge rate from the coolant tanks 
was in excess of 2500 gallons per minute. 

Instrumentation 

Two sets of thermocouple probes were installed through the walls 
of the test cylinder. Initially copper-constantan thermocouple probes 
were installed as shown in figure 20-b for the purpose of determining 
whether or not a region of two dimensional flow existed m the region 
away from the end walls. Near the mid section of the cylinder there 
appeared to be no end wall effects and the thermal field was two 
dimensional. 

The cylinder was then instrumented with a series of thermocouple 
probes to measure temperatures across the boundary layer at an l/d = 5* 
The boundary- layer probes were made from 30 gage copper-constantan 
wire with fiber sheathing left on the probes to reduce axial-conduction 
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losses along the wires. The probes were installed as shown in figure 
22, the intent being to leave the flow field as nearly undisturbed as 
possible . At large values of time the probes could affect the flow 
recirculating from the core into the boundary layer , but the very low 
core velocities and core stratification found in the numerical solu- 
tions indicated that such disturbances should be negligible. 

The thermocouple wires passed through stainless steel sheaths 
that were silver soldered into the cylinder wall. The sheaths were 
then encapsulated with rubber sealant, and the tank was vacuum leak 
tested as previously described. The thermocouple leads were housed 
in a controlled ice point cold junction box. The output wires from 
the cold junction were then lead to an analog to digital data- 
recording system where they were read out on computer tape . Care was 
taken in locating the boundary-layer thermocouple probes. These fine 
wires were quite easily displaced from a given position so that after 
being located in the correct position no further instruments or 
probes were inserted into the test cylinder. Also charging and 
purging of the cylinder was done at a slow rate to minimize convection 
velocities. Pressure measurements within the cylinder were made with 
a 0 to '15 psi absolute Statham gage and for the high Grashof number 
test with a 0 to 100 psi Statham gage. These gages were calibrated 
prior to each series of tests using a Wallace and Tiernan absolute 
gage as reference. The coolant temperature was monitored both in 
the coolant tanks and in the annular space surrounding the test 
cylinder by thermocouples that also lead through the cold junction to 
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the data recording system. 

Test Procedure 

The test procedure consisted of charging the cylinder with dry 
air at ambient temperature, sealing off the inlet and outlet valves to 
the cylinder, and then allowing a two hour waiting period for convec- 
tion currents to damp out within the cylinder. Thirty minutes prior 
to a test, water and cracked ice were introduced into the coolant 
tanks, and the coolant was brought to a uniform and steady temperature 
that typically was 1+94°R. With all thermocouples displaying steady 
state readings the plug valves were opened, and the coolant was dumped 
over the test cylinder walls. Immediately upon opening the plug 
valves , a circulating pump was activated to minimize temperature 
gradients in the cold liquid by drawing cold liquid out of the 
annular region surrounding the test cylinder, and spraying the coolant 
back over the ice crystals in the coolant tanks . Strainers at the 
bottom of each tank prevented solid ice from going past the plug 
valves and into the annular tank chamber. 

For a period of approximately five minutes following the coolant 
dumping, data sampling of all thermocouples and pressure instrumen- 
tation was taken. The digital system sampled each channel 400 times 
per second and stored the values on tape. A printer output from the 
digital system gave separate channel print outs at a rate of about 5 
channels per second for convenience in visually following the temperature 
and pressure changes from the digital system. 
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Following a test, the coolant was removed from the system and room 

temperature water was flushed through the coolant tanks, valves, and 

annu lar chamber. Dry, ambient air was pumped through the test cylinder 

and it was sealed off when all thermocouple readings showed steady , 

ambient temperatures existed within the cylinder. A two hour waiting 

period, was then used as previously described. 

Two series of tests were made using Freon and dry ice as the 

coolant . Wall temperatures down to 400°R were achieved but it was found 

almost impossible to maintain uniform coolant temperatures . The problem 

was due to the formation of "snow" when dry ice was sublimed in Freon. 

A more reliable approach for raising the Grashof number was taken by 

pressurizing the cylinder. The experimental data reported here was 

obtained at a nearly constant gas to final wall temperature ratio of 
T 

T*= .936. 

l i 

Attempts to Measure Gas Velocities 

Considerable effort was put into an attempt to measure gas velo- 
cities near a wall in natural convection flow. Early studies by 
Martini and Churchill (reference 17) using titanium dioxide dust were 
difficult and with considerable uncertainty. The attempt made here 
involved the use of helium filled soap bubbles • A flat plate was 
constructed to test against previous experiments and theory. This 
plate is shown housed in a plexiglass cage in figures 23-a and 23-b. 

At time zero hot water was forced through three passages internal to 
the plate. At the same time, bubbles that were filled with mixtures of 
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He and N 2 and were neutrally buoyant were introduced to the lower 
leading ed^e region of the plate. Color movies were taken of the rapid 
entrainment of the bubbles into the plate boundary layer, and their 
subsequent acceleration vertically upward. Time displacement studies 
were made of the bubbles. Almost independent of bubble size, the bubbles 
all seek out a single strea m line in the plate boundary layer. Thus it 
appears nearly impossible to measure a velocity profile. Ihe reasoning 
for why the bubbles behave in this manner is based on the fact that the 
boundary layer is quite thin for the plate tested, and thus a steep 
velocity gradient normal to the plate exists which in effect produces 
a gradient across the bubble. Such a velocity gradient would act to 
draw the bubble laterally inward to a position of peak velocity within 
the boundary layer. A comparison of the measured velocities for both 
transient and steady flow conditions indicated that the bubbles were in 
fact traveling with close to the maximum theoretical velocities deter- 
mined by Siegel (reference 25). The solution to the difficulty appears 
to be in producing a very thick boundary layer such that the bubble 
size is small compared to the change in velocity across a distance of 
one bubble diameter. Extensive tests were made to produce extremely 
small bubbles but below about .05 inches diameter neutral buoyancy 
cannot be achieved at standard atmospheric conditions. The hope of 
mapping out the velocxty dxstrxbutxon within the cylxnder by tracing 
the bubble displacements had to be abandoned. 



Discussion of the Experimental Results 


Chapter VII 

Comparison of the numerical results and the experiments necessitated 
a consistent formulation of "both the dimensionless variables as well as 
identical wall boundary conditions. Also the numerical solutions must 
have the initial conditions imposed at the same time zero as actually 
occurred in the experiments. It was found that within 0.4 sec after 
application of the cold liquid to the cylinder wall, temperature drops 
were detected within the gas near the wall. Thus time zero could be 
determined within 0.4 second by thermocouple signals alone. For 
comparison with the numerical results , the heat conduction equation was 
used for determining time zero for the experiments by allowing time 
zero for the conduction problem to occur when the liquid was first 
dumped from the hoppers. Time zero for the numerical computations and 
the fluid was taken when a 0.2° Rankine decrease in wall temperature 
had occurred due to conduction. (See equation 5-9) 

As discussed in chapter 5, the near singularity in for a 

time dependent wall temperature, limits the value of the constant, a, 
in equation 5*9 to values less than about .99. 

With these considerations the measured and computed temperature 
distributions for an azimuthal angle of 90° are shown in figure 24. 

These results correspond to case II of table I with the exception of 
a variable wall temperature as mentioned above. The numerical solutions 
for 1.2, 4.0, and 8.0 seconds of real time compare favorably with the 



1+2 


measured results. The large computing time for the program prohibits 
carrying out the numerical solutions to larger values of time . The 
gradual thickening of the thermal houndary-layer is seen in these 
figures. Because the core acts as a reservoir of warm fluid for supply- 
ing the boundary layer, the temperature distributions within the layer 
tend to retain their profiles over a long period of time after the 
velocity field starts decaying. To illustrate the temperature decay, 
figures 25 show typical plots of both the experimental and theoretical 
temperature at a selected radial location. Figure 25”& shows an early 
time behavior both from experimental measurements, and from the numerical 
solutions that clearly indicates the inflection produced by the flow as 
the velocity distribution shifts toward a more developed profile. 

Similar behavior is seen in figure 25-b at 0 = 70 . Finally at the 
bottom of the cylinder where the azimuthal motion of the fluid ceases, 
the temperature decay appears as shown in figure 25~c • The numerical 
results predict a more rapid decay than was actually measured. Several 
possible reasons for the disagreement are available. The most likely 
source of error lies within the framework of the differencing scheme 
near the mid-plane of symmetry. Because both the azimuthal and radial 
velocities are negative near the bottom of the cylinder, forward 
differencing is used throughout the momentum and energy equations. But 
forward differencing carries the grid points into regions of largest 
change and tends not to balance with the backward grid points that, 
for the present physical system, lie in regions of lesser change. Near 
the mid— plane the azimuthal velocity goes to zero but forward differencing 
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tends to override this due in part to the finite grid size, and thus it 
appears that errors may be largest near the bottom mid-plane of symmetry 
In this respect the solutions are dependent upon the cylindrical geo- 
metry being considered. This represents a limitation or at least an 
undesirable aspect of one sided finite difference techniques. 

It is of some interest to observe the experimentally measured 
pressure decay within the cylinder. This decay represents a three- 
dimensional (closed volume) phenomenon that is in a sense primarily 
dependent on the two-dimensional flow field that transports heat to the 
cylinder walls. Figure 2 6 shows a typical measured pressure decay for 

7 

the case of a maximum Grashof number of T x 10 . The pressure decay 
is a direct measure of the over all energy loss from the fluid. 

The mass in the enclosed cylinder remains constant over the total 
volume and thus values of the average tank pressure are a function of 
the average fluid temperature within the cylinder. 

The single most important comparison made between the numerical 
results and the experiments is that of the Nusselt-Grashof relation 
for the time dependent flow. Figure 27 shows the measured time 
dependent variation of this dimensionless grouping. The comparison is 
favorable. Many steady state natural convection flows can be described 
by the Nusselt-Grashof correlation and from figure 27 it appears that 
a time dependent correlation could be written for an internal flow 
with the boundary conditions presently under consideration. The 
formulation of such a relation was not attempted in this study . 
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The Nusselt-Grashof relation is asymptotic with zero at large 
time as seen by: 



When T •+• 00 we find 4>_ 0 *► 1 and the right hand side of equation 7-el 

J »*lim-1 

goes to zero. The value of H will pass through Helium's steady 

state value (0.326) about l4 seconds after initiation of the flow. 

If figure 27 is compared to figure 18 the effect of a steadily decaying 
wall temperature is apparent. The curve shown in figure 27 can be 
interpreted as being a result of lower values of heat transfer when the 
wall is cooled unsteadily than when a step function cold wall is applied, 
such as for the figure 18 conditions. 

Examples of typical experimental radial temperature distributions 
are shown in figure 28. In figure 28a at t = 4.0 seconds, the wall 
has cooled to about .983 of the initial temperature, and thus the 
distribution terminates at the wall (R=1.0) in the manner shown. 

The thermocouple located at R = 0.2 was considered essential for 
monitoring the core temperature. If a very low ratio of T had 
been utilized, it could be anticipated that fluid would fall downward 
from the upper most walls of the cylinder. In such a case a negative 
vertical temperature distribution within the core might have been 



observed. For the T /T. ratio of these experiments the flow was 

W 1 

almost completely confined within the wall boundary layers, and no 
negative upward distributions were measured. Figure 28-b shows the 
temperature distribution at 8.0 seconds after flow initiation. The 
measured profile has broadened more than the numerical solutions 
predict. Figure 29-a shows a similar result for 0 = 70°. 

A series of experiments at Gras ho f number values up to about 
1.8 x 10^ were made and some results are shown in figures 30. The 
very slow time development of the flow field is indicated by both the 
numerical calculations which take almost an order of magnitude longer 
computational time, than the case II solutions and by the experiments. 
The agreement between the numerical solutions and the experiments is not 
as good at low values of the Grashof number as at the higher Grashof 
numbers. The very slow development of the flow reduces the early time 
accuracy because of the small differences in temperature that must be 
measured. Figure 31 is a plot of the time history of the boundary- 
layer temperature at R = .966. The slow decay is an example of the 
lower Grashof number behavior. However, the relatively more rapid 
reduction of the gas to wall temperature ratio for the lower Grashof 
number case produces an undershoot of the Nusselt— Grashof relation. 

This is due to the differences between the velocity field development 
and the thermal field development and is shown in figure 32. An 
initial response to the near singularity at early times due to low 
values of the Grashof number causes the undershoot in the numerical 
results. Thus although the disturbance is damped by the difference 
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scheme, the accuracy of the results is severely lessened by the 
destabilizing influence of a very low early time Grashof number. Re- 
stated, the low Grashof number computations are characterized by a 
relatively rapid wall cooling which raises the wall heat transfer more 
rapidly than the one fourth power of the Grashof number. The Nusselt- 
Grashof behavior in both figure 27 and figure 32 must be viewed as 
being representative of one specific flow configuration, and also a 
very specific boundary condition. The difference between figures 27 
and 32 indicates that a single correlation equation from the numerical 
results presented here, may not be possible to achieve for the entire 
range of unsteady laminar natural convection flow within a horizontal 


cylinder. 



Conclusions 


Chapter VIII 

The present investigation has revealed the dynamic and thermal 
behavior of a confined fluid subjected to gravitational body forces. 

For a semi-infinite horizontal cylinder with uniformly cold walls the 
principle features of the flow involve a rapid development of the 
boundary layer adjacent to the walls, a tendency for fluid stratifica- 
tion in the lower regions of the cylinder, and a slow decay of the 
velocity and thermal fields with time. Experiments made in this 
investigation substantiate the thermal behavior predicted by numerical 
solutions to the quasi-compressible Navier-Stokes equations for the case 
of a time dependent wall temperature decay. Because the thermal and 
velocity fields are strongly coupled, the experimental findings imply 
that the velocity field may also be accurately described by the 
numerical results. 

The following conclusions can be drawn from this study: 

1. For very low velocity natural convection flows, windward 
finite differencing, which gives first order accuracy in time and space, 
is a suitable numerical scheme when positive and negative velocities 
occur. The large computing times required for such flow fields rules 
out more time consuming differencing schemes at the present time. 

2. The windward scheme in cylindrical coordinates was not 
extremely sensitive to grid spacing ratios except when the ratio 

— > 0.2. Above this value, numerical instabilities occurred regardless 
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of the time step used. 

3. The differencing scheme appeared to lose accuracy near the 
bottom of the cylinder where both the radial and azimuthal velocities 
are negative, and the azimuthal velocity is decelerating to zero at 

0 = 0 . 

4. The early establishment of a positive upward temperature 

gradient within the central core flow, along with the effects of fluid 

viscosity combine to produce a strong resistance to induced fluid 

motion within the core. After initiation, the flow within the 

boundary layer develops rapidly, and decays over a long period of time. 

1 3P* a 3P# 

5. The dynamic pressure gradient terms — and are 

negligible in the azimuthal and radial momentum equations for the 
conditions considered in. this investigation. 

6. The relationship between the Nusselt and Grashof numbers is 
found both numerically and experimentally for the case of unsteady 
natural convection flow within a horizontal cylinder subjected to 
uniformly cold wall boundary conditions. For the case of a time 
dependent wall boundary condition, the Nusselt -Grashof relation inter- 
cepts the steady state value of the dimensionless group about l4 
seconds after the commencement of flow down the cylinder walls. For 

the case of a constant cold wall condition from time zero, the numerical 
results show that the Nusselt-Grashof relation intercepts the steady 
state value at about 3 second after commencement of the flow. 

7. No first order vortical motion was found within the core flow. 
The absence of this motion is attributed to the boundary conditions 



which require azimuthal deceleration of the boundary layer flow near 
the bottom of the cylinder. Both the theoretical and experiment al 
models establish a mid plane of symmetry that satisfies the boundary 
conditions imposed in this study. 
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Stability Analysis 


Appendix A 

The von Neumann stability analysis deals with linearized forms 
of the governing equations. The coefficient velocities are all consi- 
dered constant in equations 3»10 through 3*13. Th e effects of amplifi- 
cation of disturbances due to non linear terms is not within the scope 
of this method, however, the von Neumann criteria does provide an 
approximate stability limit that has proved quite suitable for non- 
linear fluid flow problems. The analysis will be made for the coupled 
system composed of the azimuthal momentum equation and the energy 
equation. The differences for the cross derivatives that do not appear 
in equations k.k and 1+.7 will be included for stability considerations. 
Define the following equalities: 


B = 


IuIat 

RA8 







A-l 



Use of absolute values on U and V allows the results of the 
analysis to be applied to fluid flows with either positive or negative 
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velocity components. Substituting equations A-l into the energy 
equation gives: 

4>j + J = <J)j ^(l - G - B - 2DE - 2EC ) + <(>“^(0 + EC - ~) 

♦ * ED) * * +!Ui <E0 * ?> A ' 2 

Similarly equation 3 10 yields 

C ■ ^ t 1 - 0 - B - CH (zm * 77f) - H 
+ [° + ® + + DH) 

+ 'S.w Hfrff i W ( dh) 

- sin 6(1 - 4 >j >£ )At - G A-3 


For equation A-2 we define 


a^ = l- G- B- 2DE - 2EC 


EF 

a 2 - G + EC - r 


— B + ED 
a u = ED 


A-k 


a c = 


EF 

EC + f- 



and for equation A-3 we define: 


b-l-G-B-CH 




2DH 


”2 ’ ° + °H (ftf ) 


b 3 = B + AH 


* - cn R + M 

^ R + f 


b 5 = DH 


Substituting equations A-U into A-2 gives 


Ci ■ Vj,* * * VS-l.l + a 4&l.l * a 5+5,l»l 


and substituting equations A-5 into equation A-3 gives: 


U n+1 = 
U J,Jl 


+ b 2 U J,5--l + *3^-1, i + 


+ b 5 “S-l,l 


- sin 9(1 - ♦j} il )&T - jj 

The amplification matrices formed from the coefficients in eqs. 
and A-T have characteristic values (eigenvalues), and h^. 

For the energy equation: 
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-ikpAR -ik RA0 ikgAR ikgRA0 

h l = a l + a 2 e + a 3 e + % e + a 5 e A " 8 

For the momentum equation: 

-ikpAR -ik RA0 ik AR ik RA6 

h 2 = hi + b 2 e + t 3 e + + b 5 e A-9 

The von Neumann condition can be written as |h| <_ 1.0. The a^ and 
are all real and positive except for a 1 and b.^ which may be 
either positive or negative. As indicated by Richtmyer (reference 2 6 ) 
and Heliums (reference l6) the largest absolute values of h^ and hg 
will occur when all terms in equations A-8 and A-9 are real i.e. 
when 


^RAS = kgAR = ir 


or 


k^RA0 = kgAR = 2 it 


For the maximum real value of h^ we write: 


h- * a + a. + a, + a. + a_ 
1 max 1 2 3 4 5 


A- 10 


or substituting from equation A-k we get: 


EF 

„ . __ EF + B + 2DE + CE + — A-ll 

h = 1 - G - B - 2DE - 2CE + G + CE - r- 2 

1 max * 



56 


or 


h = 1 

1 max 


A-12 


For the minimum value we find: 


\ min " a l " a 2 " a 3 " % " & 5 


A-13 


Substituting from equations A-U gives: 


h . = 1 - G - B - 2DE - 2CE -G-CE+|“-B- 2DE - CE 

1 mm 2 




or 


h . =l-2-G-2B- UDE - UCE 

1 mm 


A-15 


For stability ve must have 


2G + 2B + UDE + UCE < 1 


A-l6 


Equation A-l6 is the stability equation governing the energy equation. 


Substituting equations A-l into equation A-l6 gives: 


2lvlAt . 2lulAt UAt _ 


UAt 


AR RA 6 


(RA0) 2 P r ^T AR 2 P r ^ 


< 1 


A-17 


or rearranging gives the maximum allowable time step for stability 


of the energy equation: 
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At < 


s]vL + £lnJ. + i + i — 

ar Rie (Rae) 2 P r >^' aF|2p r '^7 


Similarly for the momentum equation we find that: 


h 2 «x = \ + b 2 * V ^ + b 5 


Substituting equation A- 5 into equation A- 19 we get; 


h 0 = 1 - G - B 
2 max 


G + OH (7tt) + B + M + OH (;7f) 


+ DH 


Simplifying gives: 


. _ CH R CH AR CH R 

h 2 max ~ R c_^. r + ~ 

R " 2 R 2 R 2 


CH AR CH R CH R 

R + M- R + ffi-R-f 


or 


h 0 = 1 
2 max 




A -18 


A-19 


A-20 


A-21 


A-22 
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If we exclude the origin then R will never he less than AR so from 

equation A-22 we see that h„ cannot exceed unity. (If R is 

2 max 

large the term in parentheses in equation A-22 will he small). Now if 
R = AR equation A-22 gives: 


h 2 max ~ ^ 



A-23 


This leads to the expression 


At < AR 2 v€^ A-2lt 

For large Grashof number flows equation A-2U does not impose a severe 
restriction on the allowable time step. For the minimum value of hg 
we write : 

h 2 min ’ h - *2 - ”3 ' H ‘ b 5 A ' 25 

or substituting equations A-5 into equation A-25 and simplifying gives 

^max’ 1 - 20 - 2 ®- A ' 26 

For stability we write: 


2G + 2B + l*DH + HCH < 1 


A-27 
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Substituting equations A-l into equation A-27 gives: 


2 |vlAt + 2 1 U 1 At + UAt + UAt < ± 

48 846 m) 2 ^~ 


A-28 


Rearranging equation A-28 gives : 


At < A-29 

2llL + iM. + J± + L_ 

AR RA6 (raq) 2 ^~ AR 2 /T 

r r 

Thus from the von Neumann approach the most stringent stability 
requirement for the present system of equations is given by equation 
A-29. (For the case where NPr <_1.0). This criteria is at most an 
approximation and requires validation by computer experiments. The 
studies reported here have in fact verified that the von Neumann criteria 
is a good approximation and any local violation of the limit given by 
equation A-29 resulted in large scale instabilities that develop quite 
rapidly in the computation. 



Thermocouple Errors 


Appendix B 

There are four sources of temperature measurement error in 
addition to actual signal readout errors that must be investigated for 
any experimental study. These errors can be classified as being due to: 
(l) fluid velocity past the sensor; (2) conduction losses along the 
thermocouple wires ; (3) radiation of energy from the thermocouple to its 
surroundings; (4) transient error due to a finite time lag in response 
caused by the thermal capacity of the sensor. These errors will be 
evaluated for operating conditions that are typical of the experiments 
carried out. Moffat (reference 27) gives a detailed discussion of the 
pertinent equations for such an analysis. 


Thermocouple velocity error 

The difference between the gas total temperature T T and the 


temperature, T_, given by the Junction is: 
J 


T m - T_ 
T J 


= (1 - 


u; 

0t r ) 2gJC 


B-l 


where : 

t j - T . 

a = Recovery Factor = „ - T 

T = fluid static temperature 
00 

U = fluid velocity 
00 

The recovery factor for wires perpendicular to the flow has been 
extensively studied and a summary of the results is shown in reference 


6 0 
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27. For the very low fluid velocities encountered in the present work 
the recovery factor can he conservatively taken as: = 0.6l. For a 

maximum fluid velocity of 4 feet/second the velocity error is 
calculated as : 


(.39)16 
2 ( 32 : 17TT1 


The velocity error appears to he negligible for the present experimental 
conditions. 

Conduction Error 

The losses due to ml conduction along the thermocouple wires 
are given hy the following equation: 



B-2 


where : 

T t = Fluid Total Temperature 
T m = Mount Temperature 

L = Distance from the Wall to the Sensor Junction 
d = Wire Diameter 

h p - Convective Heat Transfer Coefficient 

k c = Thermal Conductivity of Sensor 

s 

v = Thermal Conductivity of the Gas Evaluated at Stagnation 
G 

Conditions 

The Reynolds number based on wire diameter is given as; 
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R = ^ = = 22 .8 B-3 

ed y T 1.46 x 10~ k 

If the correlation equation of reference 27 for wires and perpendicular 
to the flow is used, the following results are obtained: 


N = (.44 + .C^R 1 ^ 2 B-4 

u “ ed 

k N 

with h = ^ 

c a 

so that ®-6 

For a fluid total temperature T^, = 536°R and a mount temperature 
= i+92°R equation B-2 gives: 

T T " T J “ cosh 4l4 


Thus 


T t - Tj « 1°R 


B-8 


The conduction losses are seen to be entirely negligible even for the 
shortest thermocouple sensors used in this investigation. 


Radiation Losses 

The radiation loss for a thermocouple within an enclosure is given by: 



k a e 
r 

h 

c 




B-9 
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where 


K = Radiation Form Factor 
r 

a = Stefan Boltzman Constant 
A 

= Ratio of Area Available for Radiative Loss to Area Available 
c for Convective Loss. 

e = Emissivity Compared to a Black body 

For an enclosure that is large compared to the wire dimensions it is 

known that K ~ 1.0. 

r ^r 

For the present problem we may also write = 1.0. Also for oxidized 


c 

base metal thermocouples e * .85. And 


= 4.75 x 10' 




Ft 2 sec°R 


B-10 


Considering the largest gas to wall temperature difference to be 20°R 
for the unsteady wall temperature decay of the present experiments, 
equation B-8 gives: 


T - Tj = 3.70 jc 10 -11 (1.1U ic 10 10 ) 

Thus (Tj ' ’ °- k2 ° R 6-11 

Thus the radiative losses for the case of an unsteady wall temperature 
decay can be of significance. Within the boundary layer there is 
some reduction of the loss due to lower fluid temperatures. This 
affect is partially overcome by a reduction in the value of h^ within 
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■boundary layer. As a percentage of the total gas temperature, equation 
B-ll represents a very small correction. It is however a known source 
of error and equation B=8 must be applied to the thermocouples that 
reach a junction to wall temperature difference near the values used to 
give equation B-ll. 


Transient Error 

The lag in thermocouple response can be characterized by: 




B-12 


where 



B-13 


w is a characteristic time that depends directly on the thermal 
storage capacity of the thermocouple and inversely upon the heat input 
to the thermocouple per degree of temperature difference. 

For Copper 


C = .0918 


BTU 
Lb °R 


so that 


b> = i2i. = .751* seconds B-lU 

4(1.09 x 10" 2 ) 12 


Thus 
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= 0.55 °R/sec . 

\dt / max 

The maximum transient error is close to the radiative error and very 
neatly compensates for the radiative error calculated previously. 

Other errors were reduced to a minimum for the system by calibra- 
tion of the pressure instrumentation prior to each series of tests 
and also by performing a precise balance of the readout system prior 
to each test. 

A large number of tests were carried out and each test condition 
reported here was duplicated at least three individual times. The 
comparison of temperature and pressure histories for duplicate tests 
is very close with random variations in temperature readings of no 
more than 1 °R between two tests run at the same conditions. The time 
history of pressure and temperature within the horizontal cylinder 
could be duplicated within 0.5 °/o from one test to another. This 
close repeatability was obtained over the entire Grashof number range of 
the experiments. Systematic errors are established by the limits of 
output accuracy of the sensors and pressure transducers and the total 
readout system. The rapid sampling rate of the automatic readout system 
removes the need for analyzing the experimental uncertainty by means 
of a Gaussian error curve. Instead, it is seen that the fluctuations 
shrink the distribution curve inside the smallest scale divisions of 
the readout system. As an illustration data samples are taken at a rate 
of U00 samples per second. The peak rate of temperature change measured 
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was close to 1°R per second so that effectively 100 sample readings 
were taken within each 1/4 °R change in the value of t emper attire . For 
machine plotted data this is about the limit of scale size that can 
usefully he worked with. The data readings are truly scale limited 
rather than being perturbation limited. Hie absolute uncertainty of 
the temperature measurements is set by the error limits of the 
thermocouple junction itself. All known system errors are far below 
the thermocouple error which commercially is given as + 1 1/2°F. 
Calibrations of the thermocouples within a controlled ice point instru- 
ment showed less than 1/2 °F discrepancy between recorded signal and 
the ice point temperature. The temperature measurements have an overall 
uncertainty of something on the order of 1 1/2°F which is within about 
+0.2 o/ o of the measured values. 



Table I 


Parameter Variations for Numerical Studies 


Case 

I 

Case II 

Case III 

m = 

I. 

8 . 3 U x 10 

NG ■ 7 x 10 7 

NG * 1.36 x 10^ 

r 

r 

r 

II 

.715 

P = .715 
r 

P r = .715 

T 


T 

T 

II 

.936 

=£ ■ .93 6 
,T^ 

“ - -936 
i 

II 

> 

EH 

constant 

T = constant 
w 

T = constant 
w 

r v = 

0.1 ft. 

r w = .25 ft. 

r = .25 ft. 

V 

P. = 
■ 1 

1 atm. 

P i = 6 . 7 ^ atm. 

P^ = 1 atm. 


Case II - T 


Case III - T 


Identical to Case II 

except : 

Identical to Case III except: 


T = <Ko,t)(T. - T ) + T 
w 1 x o X o 


T = $(o,t)(T. - T ) + T 

V 1 * 
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Figure 3a.- Flow diagram for finite difference calculation 
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Figure 3-1 • - (Continued). 
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Figure 3-c.- (Continued) 
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Write tape 2 T, At, u, v, v? , 
Print final T, At 


m 

STOP 2 


Figure 3-d.- (Concluded). 







Figure k&.- Azimuthal velocity variation with radial distance. 







Figure U-b.- (Continued) 












Figure U-c.- (Concluded) 












Figure 58'*“ Radial velocity variation, with radial distance 















Figure 5^.- (Continued) 
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Figure 6b.- Temperature variation with radial distance 
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Figure 7.- Wall heat transfer decay with time. 














Figure 8.- Azimuthal temperature distribution within the core. 













Figure 10b.- Radial grid refinement effects at a Grashof number 1 
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Figure lid.- Nusselt-Grashof relation for lov Grashof numbers natural convection flow within a 
horizontal cylinder. 



Figure 12.- Azimuthal velocity distribution in horizontal cylinder 












Figure Ik.- Azimuthal velocity distribution for time dependent wall temperature. 




Figure 15.- Time development of the Grashof number for unsteady wall temperature. 







Figure l6a.- Temperature distributions for a time dependent wall temperature . 










Figure 17a.- Azimuthal velocity distributions for constant and variable wall 
temperature . 





Figure 17b.- Azimuthal velocity distribution for constant and variable vail 
temperatures . 
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Figure 18.- Nusselt-Grashof correlation for time dependent flow within a horizontal cylinder. 
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Figure 20a.-* Instrumented natural convection chamber. 



Figure 20b.- Schematic diagram of 



erimental natural convection apparatus. 
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Figure 21.- Cooling manifolds for natural convection chamber 
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Figure 23b*- Vertical flat plate apparatus concluded 
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Figure 24a.- Temperature distribution for a time dependent wall boundary condition. 
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Figure 24c.- Temperature distribution for a time dependent wall boundary condition. 
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Figure 25a.- Boundary layer temperature decay at 0 - 90° and R - .966. 





Figure 25c.- Boundary layer temperature decay at 0=0° and R - .966. 





Figure 26.- Measured pressure decay in natural convection flow within a horizontal cylinder. 
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Figure 27.- Nusselt-Grashof relation for natural convection flow within a horizontal cylinder 
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Figure 28b.- Radial temperature distribution at 0 - 90 , t 
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Figure 30a.- Radial temperature distribution at 0-9O,t-2.O seconds. 
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Boundary layer temperature decay at 6 = 90 and R = . 966 , G 


= 1.3 x 10 


Figure 31 
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Figure 32.- Nusselt Grashof relation for 6 = 90°, G^ max - 6.5 x 10" 
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